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A reactive fluid dissolving the surface of a uniform fracture will trigger an instability in the dissolu- 
tion front, leading to spontaneous formation of pronounced well-spaced channels in the surrounding 
rock matrix. Although the underlying mechanism is similar to the wormhole instability in porous 
rocks there are significant differences in the physics, due to the absence of a steadily propagating 
£N) ' reaction front. In previous work we have described the geophysical implications of this instability 

in regard to the formation of long conduits in soluble rocks. Here we describe a more general linear 
stability analysis, including axial diffusion, transport limited dissolution, non-linear kinetics, and a 
04 | finite length system. 
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I. INTRODUCTION 



Fracture dissolution is an important component of a number of geological processes, including the early stages 
of karstification diagenesis |2|, and the evolution of carbonate aquifers Q. It also plays an important role in 
geoengineering applications such as dam stability [|| , oil reservoir stimulation methods @ and leakage of sequestered 
■ CO2 [6| . The dynamics of evolving fractures is complex, due to the highly nonlinear couplings between morphology, 
flow and dissolution. Theoretical p], 0, H| and experimental studies pHll| have shown that the positive feedback 
between fluid transport and mineral dissolution leads to an instability in an initially uniform reaction front and the 
subsequent formation of pronounced dissolution channels, deeply etched into the rock surfaces. These processes were 
shown to be important in the development of limestone caves |8[, and also in the assessment of subsidence hazards, 
since they dramatically speed up the growth of long conduits. Understanding spontaneous flow focusing during 
fracture dissolution is also important to the petroleum industry, for efficient acidization of natural fractures and for 
acid fracturing of porous rocks. In the former process, acid is pumped into the fractured reservoir to dissolve material 
blocking the pathways between the wellbore and the reservoir. Spontaneous channeling increases the effectiveness of 
the process by creating highly permeable pathways, minimizing the amount of acid needed. In acid fracturing the 
fluid pressure is high enough to induce hydrofracturing; the newly created fractures are then etched with acid to 
increase the permeability of the system. Nonuniform dissolution is crucial in this process, since a uniformly etched 
fracture will close tightly under the overburden once the fluid pressure is removed; significant permeability will only 
be created by inhomogeneous etching when the less dissolved regions act as supports to keep more dissolved regions 
Q\ ; open. 

In this paper we investigate the initiation of the instability in a fracture dissolution front and assess the wavelength 
and growth rate of the most unstable mode as a function of physical parameters characterizing the rates of transport 
and reaction in the fracture. In Sec. |TI]we present the two-dimensional averaged equations for fracture dissolution; a 
^v^j . detailed justification of the transport equation (|3]) is given in Appendix Next we consider a uniform fracture where 
an analytic solution is possible; this forms the base state for the subsequent stability analysis in Sec. IIVI Results are 
presented in Sec. IVI1 extending our previous analysis 8] in several directions. We now consider axial diffusion of 
reactant as well as lateral diffusion and also the effect of cross-aperture diffusion on the effective reaction rate. After 
that we lift the assumptions that the fracture is of infinite length and that the reaction kinetics are linear. We finish 
with a summary of our results and conclusions. In a subsequent paper we will describe an analysis of the instability 
in the dissolution of a porous matrix. 



II. EQUATIONS FOR FRACTURE DISSOLUTION 

Fractures are geometrically characterized by a short dimension (z direction), the aperture, and two much longer 
dimensions, length (x direction) and width (y direction). In natural fractures the aperture is typically less than 
1mm, while the length (L) and width (W) are of the order of meters (see Fig. [T]). It is typical to exploit this 
difference in scales by introducing approximate two-dimensional equations for fluid flow, reactant transport, and 
erosion. Fluid flow is described by the Reynolds equation for the local volume flux (per unit length across the 
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FIG. 1: Dissolution of a rough fracture of length L and width W; fluid flow is in the x direction and the fracture surfaces 
dissolve in the normal (z) direction. The aperture h(x,y) is the distance between the fracture surfaces. 



fracture), q(x,y,t) — L v(x,y,z,t)dz: 



h 3 

q = - — V Pl V-q = 0, 



(1) 



where /i is the fluid viscosity. The essence of the Reynolds approximation is to assume that the exact result for 
stationary flow between parallel plates can be applied locally to a varying aperture. In this approximation the 
pressure is independent of height and reduces to the two-dimensional field p(x,y). The validity of the Reynolds 
approximation for rough fractures has been examined in [l2j and [l3j]. The key requirements are: (i) low Reynolds 
number flow, Re <C 1 (ii) slow variation in aperture |Vft| <C 1. We will assume these conditions hold in what 
follows. The incompressibility condition in Eq. ([1]) ignores effects of the reactant (or product) concentration on the 
mass density of the fluid. This assumption is valid for the majority of natural systems; for example, in limestone 
dissolution the density correction due to the dissolved species is of the order of 0.01%. However, dissolution of halite 
(rock salt) is a notable exception; here the increase in mass density can be as large as 25%. 

The transport of reactant can be described in terms of a two-dimensional concentration field that has been averaged 
over the aperture. The most important average is the "cup-mixing" or velocity-averaged concentration [l4j |. 



c(x,y,t) = 



1 



q{x,y,t) 



h(x,y,t) 



\v(x,y,z,t)\c 3d (x,y,z,t)dz, 



(2) 



where we use c^d to identify the three-dimensional concentration field. Under certain conditions, discussed in Ap- 
pendix [AJ the three-dimensional convection-diffusion equation for reactant transport in the fracture can be reduced 
to a two-dimensional convection-diffusion- reaction equation for the cup-mixing concentration [3, 0, HI , 



q ■ Vc = DVh ■ Vc~2R{c) 



(3) 



where R(c) accounts for reactant transfer at each of the fracture surfaces. The slow dissolution of the rock surfaces 
allows the time-dependence in Eq. ^ to be neglected (Appendix IA II) . 

In this paper we will usually assume a first-order dissolution reaction at the fracture surfaces R = kc w , where k is 
the rate constant and c w is the reactant concentration at the fracture surface. The reactive flux R must balance the 
diffusive flux at the surface 



Rdiff = -D{Vc) v 



(4) 



where the gradient is pointing towards the surface. Alternatively, and more usefully, the diffusive flux can be expressed 
in terms of the difference between the surface concentration, c w , and the cup- mixing concentration, c by using a mass- 
transfer coefficient or Sherwood number }l4| . 



Rdiff 



DSh 
2h 



(c- Cw). 



(5) 
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The Sherwood number, Sh, depends on reaction rate at the fracture surfaces (k) but the variation is relatively 
small [TBI, [ill], bounded by two asymptotic limits: high reaction rates (transport limit), Sh = 7.54, and low reaction 
rates (reaction limit), Sh = 8.24. In the numerical calculations we approximate the Sherwood number by a constant 
value Sh = 8. 

By equating the reactive and diffusive fluxes R = Rdi //we obtain the standard relationship between c w and c , 

(6) 



1 + 2kh/DSh 

The reactive flux can then be expressed in terms of the cup-mixing concentration, 

R ( C ) = k effC, (7) 



where the effective reaction rate is given by 



keff{h) = i + 2fchAPSh ' (8) 



In sufficiently narrow apertures the dissolution kinetics are reaction limited and the concentration field is almost 
uniform across the aperture so that k e ft ~ k. However, as the fracture opens the reaction rate becomes hindered 
by diffusive transport of reactant across the aperture. When kh/ DSh 3> 1, dissolution can become entirely diffusion 
limited with k e ff ~ DSh/2h. 

A derivation of Eq. ([3]), with the kinetics described in Eqs. ([7} and ©, will be given in Appendix [Al starting 
from the full three-dimensional transport equations. In particular, the diffusive term in Eq. ([3]) is shown to be purely 
molecular for either convective (q/D — > 00) or reaction-limited (2kh/D — > 0) transport. In taking the Sherwood 
number to be independent of the distance from the inlet, we are assuming that entrance effects are negligible. For a 
flat plat geometry the entrance length scale Z,- n is given by [13] 

tin = 0.016^, (9) 

taking U n as the distance over which the Sherwood number is within 5% of its asymptotic value. This length is small 
compared to the reactant penetration length under the typical conditions of fracture dissolution (see Sec. IIII[) . 

Equations ([7]) and <JSj> describe a dissolution reaction controlled by the concentration of reactant; a typical example 
is dissolution of fractures (or porous rocks) by a strong acid. However, when calcite is dissolved by aqueous CO2 at 
pH values similar to those of natural groundwater, the dissolution rate is limited by the calcium ion undersaturation 

Csat — < 



R(Cca) = -k e ff(c S at - C ca ), (10) 

where c ca is the flow-averaged concentration of dissolved calcium ions. The sign of R accounts for a dissolution flux 
into the fluid rather than a reactive flux into the surface and so the transport equation for the undersaturation takes 
the same form as ([3]). In the rest of the paper we will use c to represent either the concentration of reactant or the 
undersaturation of dissolved minerals. 

A reactive fluid with an inlet (x — 0) concentration ci n dissolves the surrounding rock, increasing the fracture 
aperture at a rate 

d t h = 2k eff y—, (II) 

where 7 = c in /vc so i is the acid capacity number or volume of solid dissolved by a unit volume of reactant. Here 
c so i is the molar concentration of soluble material and v accounts for the stoichiometry of the reaction. Mineral 
concentrations in the solid phase, are typically much higher than reactant concentrations in the aqueous phase and 
the characteristic dissolution time, 

t d = h/2k efn , (12) 

is large for natural minerals in typical groundwater conditions; for limestone fractures it is approximately 2 months 
Q . Thus there is a significant separation between the dissolution time scale and the relaxation of the concentration 
field (t ~ h 2 /D), which justifies dropping the time dependence in Eq. (JU); for further discussion see Appendix IA II 
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III. CONCENTRATION PROFILE IN A UNIFORM FRACTURE 

Let us first consider a uniform aperture h(x, y) = ho and find the corresponding concentration profile; the solutions 
will form the base state for the stability analysis. The flow rate qq is independent of space and the transport equation 
is 

qo d x c-Dh d?c=- T ^, (13) 
where we have absorbed the transport correction into a single factor, 

For an inlet concentration Ci n , Eq. (1131) has an exponentially decaying solution, 

c(x) = c in e~ RX , (15) 

with a penetration length l p — k" 1 given by 



The Peclet number, 

Pe=|, (17) 

measures the relative magnitude of convective and diffusive transport of solute, and the effective Damkohlcr number, 

2k eff h 2kh 

Da e// = — = 7, , r , . ( 18 ) 

qo (i- + u)qo 

relates the effective surface reaction rate, Eqs. ([7]) and ©, to the rate of convective transport. 

It will be convenient to frame our results in terms of the transport correction G (1141) and the convective parameter 



H = ^ll. (19) 

A discussion of the natural length scales of the problem and their relation to H can be found in Appendix |Bj The 
inverse penetration length can be written in terms of H, 

Pe 



nh Q = — (yi + 4ff - 1 J , (20) 



/i/i = Da e// , (21) 



with the important limiting cases: 

• convection dominated (H — > 0) 

• diffusion dominated (H — > oo) 



In Appendix E] we show that (J5J) is valid for all G when H = (Sec. IA 21) and for all H when G < 1 (Sec. IA3j) . 

For long fractures, the reactant penetration length is the natural length scale for dissolution. On the scale of k^ 1 
the entrance length is 



GSh 



nh = VPeDa e// = A / T -^, (22) 



Kl in = 0.008Pe 2 (v / l + 4 J ff- 1). (23) 

In the convective (H — > 0) limit, Kli n — 0.016GSh/(l+G) < 0.12 over the whole range of reaction rates; it is vanishingly 
small in the reaction (G -t 0) limit. In the diffusive (H -t oo) limit nl in = 0.016Pe^GSh/(l + G) < 0.05Pe, which 
is again small (since Pe <C 1). In Sec. IVIEl we will examine the instability in finite-length fractures kL < 1, but only 
in the reaction limit (G — > 0), in which case h n /L 0, even for finite kL. 
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IV. LINEAR STABILITY ANALYSIS OF A UNIFORM PROFILE 

The discussion in Sec. [HI supported by the derivations in Appendix [XJ leads to the following average equations for 
the concentration, aperture and flow fields in an evolving fracture: 

«,a,c + q ,a,c - d mhd.c) + d s (hd s c)} = ~d,h 

2k^c 

Cmdth = 1 + 2kh/DSh (er0Si0n) (25) 

d x q x + d y q y = (continuity) (26) 
3 3 

d y q x - j^qxdyh = d x q y - -q y d x h (compatibility) (27) 

Here the Reynolds equation (JTJ) has been replaced by the more convenient equations for continuity (|26l) and compat- 
ibility (|27[) (see Appendix [Cl) . When supplemented by appropriate boundary conditions: 

c(x = 0,y,t) = Ci n , c{x -> oo, y,t) = 0, (28) 

q x (x -> <x,y,t) = q Q , q v (x = 0,y,t) = 0, (29) 

Eqs. (j2"4")l ([2~T| form a complete, albeit approximate, description of the erosion of a single fracture (in the domain 
x > 0). The constant pressure condition at the inlet has been replaced by the boundary condition q y (x = 0) = 0. 

The above equations allow one-dimensional solutions in which the fields depend only on x and t. This corresponds 
to uniform dissolution of the fracture, an assumption still commonly found in models of fracture dissolution H3] ■ 
For example, in the reaction- limited, convection-dominated case (G — > 0, H — > 0), the solution is 

c(x,t) = c m e- 2kx l\ (30) 
h(x, t) = h + 2kjte- 2kx/q , (31) 

q(x,t) = q e x . (32) 

In H we showed that the solution represented by Eqs. (|30|) -([32 |) is unstable to infinitesimal perturbations along the 
y direction. Here we will not limit ourselves to the reaction-limited, convection dominated regime, but consider more 
general kinetics and transport. Thus K will no longer be equal to 2k /q, as in (|3H|) and (|3"Tj) . but instead it will be 
given by the general expression f|20[) . 

An important detail in the stability analysis is that the base state for the aperture (|3ip is itself time-dependent. 
The stability of nonautonomous systems is in general a difficult problem [2l[ and in [f| we adopted an approximate 
approach [22j in which the base state is frozen at a specific time, to, and the growth rate is then determined as if the 
base state were time- independent (the quasi-steady-state approximation). The validity of this approach was tested 
by comparing the results of the quasi-steady-state approximation with a numerical solution of the complete system 
of equations (fM)) - (|2"7| . In particular, we were able to show that the most relevant instability is obtained by freezing 
the base state at to — and in the present paper we will focus on this case. The solution at t = is 

c b (x) = c m e~ KX , h b (x) = h , q b (x) = q e x , (33) 

which simplifies the subsequent calculations. 

The linear stability analysis proceeds by considering infinitesimal perturbations to the base profile (|33p : h = h b + 5h, 
c = c b + Sc and q = q b + Sq. This gives the following linearized equations for the aperture, concentration and flow 
fields: 

SqAc b + q b d x Sc - D [h b d 2 Jc + h b d 2 Sc + 6hd 2 x c b + (d x 8h)(d x c b )] = -^d t 6h, (34) 

, 2kh A aX u , (, , 2fc/ tb y 1 {2kf 1Cb xu 
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d x 8q x + d y Sq y = 0, (36) 



3 

d y Sq x - —q b d y Sh = d x 5q y . (37) 
rib 



Terms in d x hb have been omitted from Eqs. (|34p and (I37p . since the expansion is about an x— independent aperture 
field. In Eq. (|35|) we have made use of the erosion equation for the base field, Cj„(l + 2khb / D Sh)d t hb = 2k^Cb- 

The linearized equations for fracture dissolution can be simplified by transforming to dimensionless variables. We 
take the penetration length k^ 1 as the unit of length, and the characteristic inlet dissolution time, td (|T2|) . as the unit 
of time. The dimensionless variables are then: 

* = v = Ky, r = JT — w -. (38) 

The concentration is scaled by the inlet concentration a n , while the aperture and flow rate are scaled by their 
(constant) values in the base state: 

c= — , h=—, q=— . (39) 

Cm iIq qo 

The dimensionless base-state solution is: 

c b = e~Z, h b = 1, q b = e ? , (40) 
and the dimensionless perturbations can be found from the following equations: 

^-—d T 8h = e-*5q* - d ( 6c + (d?8c + dl5c + e^Sti - e^d^Sh) , (41) 

d T Sh + -^—e-^Sh = 5c, (42) 
1 + G 

Sf <% + tfSqt = 3fl$}<&. (43) 

In deriving (I43[) we have combined the continuity equation (|36p and the compatibility equation (1371) to eliminate Sq^. 
The transport equation (|4ip involves two new dimensionless constants, each one based on the penetration length 

Pe - qa - — f44) 

2k ef f 2H , , 

Da K = — S±£ = - . 45 

q n VT+4H-1 

Pe K is the ratio of convective to diffusive fluxes on the length scale while Da K is the ratio of convective to reactive 
fluxes on the same scale. The physical significance of these parameters is discussed in Appendix |B] Rewriting the 
transport equation in terms of Pe K and Da K and rearranging to isolate the term in 8q^ , 

Sqt = e« [Da K <9 r + Pe" 1 ^"*] Sh + e« [<% - Pe^df + flg)] St. (46) 

Assuming that the perturbations are sinusoidal in r\ and exponential in r, 

8c = / c (0co B (u» 7 )e* T J (47) 
8h = f h ($cos(uri)e 0T , (48) 
Sq* = f q (Ocos(u V )e^. (49) 

Note that Cj and w are dimensionless quantities related to the instability growth rate u> and wavelength A by the 
relations 

u) = ujt d , u = — . (50) 



7 



Substituting the expansions (|47|) - (|49|) into Eqs. (|46|) . (1421) . and (|43|) leads to coupled equations for the one- 
dimensional fields / c (£), and /</(£) : 



/, = e« [Da K (i + Pe^e^] + e« [d 5 - Pe" 1 ^ 2 - u 2 )] / c 



(51) 



Ge-g 
1 + G 



/ft — /c 



(52) 



(5| - w 2 )/ 9 = -3u 2 f h . 
Eliminating f c , we express / q in terms of only 

/, = £ ( [Da K £ + Pe-^e-^] + - Pe" 1 ^ - u 2 )] 



1 + G 



fh • 



and, substituting into (|53[) . obtain a fourth-order equation for the £ dependence of the aperture field, 



(a| - zi V [Da K d; + Pe-^ce-^] + [3 4 - Pe^df - u 2 )] 



1 + G 



//» + 3u 2 / ;i = 0. 



(53) 



(54) 



(55) 



The boundary conditions on the perturbations can be found from Eqs. (|28|) and (|29|) . From the inlet and outlet 
conditions (|2"8j) it follows that dissolution at the inlet is uniform (because c = 1), 



A(f = o) = o, 

and that far downstream the aperture is unperturbed, 

/ ft (£^oo)=0. 



(56) 



(57) 



The boundary conditions on the flow (|29[) also impose conditions on fh through Eq. (1541) . The uniform pressure at 
the inlet leads to a condition on q^, 



/ 9 (£ = o) = [0*/ 9 ] £=o = o, 

which, by means of (|54j) . imposes a third-order boundary condition on fh, 
d e efi ( [Da K w + Pe" 1 ^ - *] + [<% - Pe" 1 ^ - u 2 )] 



Ge-Z 
1 + G 



= 0, 



4=o 



The outlet condition 



/,(e-^oo) = o, 



imposes a further restriction on fh, through Eq. (|54p . namely that it must decay at least as fast as e ^, 

e c /h(£ -> oo) = A 



(58) 



(59) 



(60) 



(61) 



In most cases the constant A must be zero in order for (f60| to be satisfied, but in the convective limit (H = 0), 
the solution fh = Ae~^ is an eigensolution of (j54"]) with zero eigenvalue, and therefore satisfies the far-field boundary 
condition on f q . 

Since the initial amplitude of the instability is arbitrary, the four boundary conditions impose an additional con- 
straint which can be used to solve for the eigenvalue G)(u). We have used a spectral method, which we summarize in 
Sec. El to find the dispersion relation numerically. In certain limiting cases further analysis is feasible; we describe 
these on a case by case basis in Sec. I VII 
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V. SPECTRAL METHOD 



The solution of equation ([55]). together with the boundary conditions (|56|) , (|59l) . and (jBlj) . was obtained using the 
pseudospectral, boundary-bordering method [23l. [2^|. For a given linear operator, 'H, the differential equation 

nf(0 = g(0, o<e<cx., (62) 

is represented as a linear system 

Hf = g (63) 
where the elements of the vector / are the coefficients of the expansion of /(£) in the basis functions 

N 

/(0 = £/^-i(0- (64) 

7=1 

Matrix elements of H are calculated at N — 2 collocation points, 

tf*faj = [«*i-i(Ok=* ( 65 ) 
and the corresponding elements of the right-hand-side vector are 

gi+2 = g{ii)- (66) 

The first two rows of H are used impose the boundary conditions at £ = 0. If the boundary conditions are expressed 
in terms of the linear operators Bi< , 

B v {f) = a v , i' = l,2, (67) 

then in the matrix representation 

H v j = [Bi'*i-i(0] e=0 . ft' = ( 68 ) 

where £' = 1, 2. 

The basis functions are rational Chebyshev functions in TZ + = [0, oo], defined as 

where T n (t), with n = 0, 1, 2, . . ., is a Chebyshev polynomial of the first kind, defined in the range —I < t < 1. The 
convergence of the solution depends on a suitable choice of the mapping parameter, L, which varies somewhat with 
wavelength. For small numbers of basis functions (N < 20), we took L — 1 at short wavelengths (u > 1) and L = 10 
at long wavelengths (u < 1). However, for larger numbers of basis functions (N > 50), a constant L — 10 was suitable 
for the whole range of wavelengths, 0.01 < u < 10. For a given L and N, the N — 2 collocation points are [23| . 

& = £cot2 (l7f^)> i = i,...,J\r-2. (70) 

The dispersion relation can be found by solving the linear system of equations represented by (I65|) - (|68p . with 
boundary conditions /(£ = 0) = ([56]) and = 0) = 1, which fixes the amplitude of the perturbation. Then, 

we iteratively seek the largest value of Cj for which the boundary condition in (IBTJl) is satisfied and hence find the 
dispersion relation u>(u). There is no need to separately impose the far-field regularity conditions, Eqs. (|57l) and (|61l) . 
since this is automatically incorporated by the basis functions (23[. We have cross-checked the spectral code with 
analytic solutions in a number of special cases (see Sec. IVI|) . and a Maple version of the spectral code is included in 
the Supplementary Material. 



VI. RESULTS 



In general, the dispersion relation ()55[) must be solved numerically; for example, using the spectral method described 
in Sec. [V] However, in the important limiting case of convection-dominated (H — > 0), reaction-limited (G — > 0) 
dissolution, it is possible to obtain a tractable analytic dispersion relation, as shown in Sec. I VI Al We can also obtain 
analytic solutions in other limiting cases, but the solutions are too lengthy to be reproduced in print, although we 
include Maple workbooks as Supplementary Material. Analytic calculations from Maple [25] and Mathematica (26j 
were crosschecked with each other and with the spectral code (Sec. |VJ in many cases. 
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FIG. 2: Growth rates of the inlet instability in the purely convective case (H — 0). The solid line corresponds to the reaction- 
limited case (G — 0), whereas the dash-dotted curve corresponds to the diffusive limit (G = oo) and the dashed curve is for 
mixed kinetics (G = 1). The dimensionless growth rate w = utn, Eq. (|12[) . is plotted against the dimensionless wavevector, 
u = 2-7T / k\. 



A. Convection-dominated dissolution: H 0. 

In convection-dominated flows (H — > 0), the Damkohler number on the scale of the penetration length Da K = 1, 
and the corresponding Peclet number Pe K — > oo. The dispersion relation (|55[) then simplifies to 

}A + 3u 2 f h = 0. (71) 

There is an analytic solution of Eq. (1711) in terms of a linear combination of three generalized hypergeometric 
functions z a [z — l)3i r 2({ai J 0-2, 0,3}, {61, &2jS 2), where and 6fc are complicated algebraic functions of G and m, 
z = —Go) -1 exp(— £)/(l + G), and a is a simple function of u. As the solution is lengthy and not very informative we 
do not include it here, but a Maple notebook is included as Supplementary Material. 
A much simpler equation is obtained in the reaction limit (G — ► 0) of Q71[) [8], 

(9| - u 2 )6j e Z(l + d 6 )f h + Zu 2 f h = 0. (72) 

The general solution of ([721 is 

/ h (0 = Ae~* F 2 (1 + u, 1 - n; to-^e"*) 

+ Be ( ™- 1)? 0-F2 (1 + u, 1 - 2u; SuP^e - *) 

+ Ge-( fi+1 ^ F 2 (1 + u, 1 + 2u; aST^e - *) , (73) 

where A, f?, and G are constants and qF%(j>, q;z) is a generalized hypergeometric function. The far field boundary 
condition (|6ip requires that B = 0, while the condition /^(0) = ((SBJ) is then sufficient to determine the function 
fh(£) to within an arbitrary constant, which is the initial amplitude of the perturbation. Imposing the final boundary 



(91 



1 + G 
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condition (I59[) gives a dispersion relation for Cj(u), 



(b z Q F 2 (1 + u, l + 2u;3& L u 



3(1 + 2u)Cj F 2 (2 + u, 2 + 2u; 3t2; _1 w 2 ) 



+9m 2 0^2 (3 + u, 3 + 2u; Sw^u 2 ) 



F 2 (l + u,l-u;3u]- 1 u 2 ) = 



! F 2 (2 + u,2-u;3lj~ 1 u 2 ) + 

3u 2 F 2 (3 + u,3-u;36j~ 1 u 2 ) F 2 (l + u, 1 + 2ii; Zuj^u 2 ) , (74) 

where oF 2 (p, q; z) — oF 2 (p,q; z)/T(p)T(q) is a regularized hypergeometric function |27j |. The maximum growth rate 
(largest positive root) at each u from ([71)1 corresponds to the solid line (G = 0) in Fig. [5] The positive growth 
rates show that the front is unstable across the whole spectrum of wavelengths, with a well-defined maximal growth 
rate, 0J ma x = 0.7947 , at a wavelength X max — 4.74k -1 . An individual fracture will therefore develop a strongly 
heterogeneous permeability during dissolution, with an inherent length scale that depends on the kinetics and flow 
rate (via k), but not the initial topography. There is no lower limit to the reaction rate for unstable dissolution if the 
scale of the fracture is sufficiently large. 

Figure also shows the impact of reaction kinetics (controlled by the parameter G) on the dispersion relation. 
For wider apertures (i.e. C > 1), diffusional transport of reactant across the aperture has a stabilizing effect on 
the growth of the instability. The fastest-growing wavelength, X ma x, is pushed towards longer wavelengths and at 
sufficiently short wavelengths perturbations in the front are stable. 



B. Reaction-limited dissolution: G — > 0. 



The dispersion relation (1551) can also be solved analytically in the reaction limit, G — 0; the solution of the dispersion 
equation, 

(d 2 - M 2 )e« { [Da«w + Pe~%e-t] + <2> - Pe^df - u 2 )] } f h + 3u 2 f h = 0, (75) 

is again a combination of hypergeometric functions 3^3 ({oi, a 2 , 03}, {b\, b 2 , 63}; z) exp(g£), where et^, bk and g are 
functions of u and H, and z = — exp(— £). Again we have included a Maple notebook in the Supplementary 
Material. 

In the diffusive limit (H — > 00), Da K — > Pe -1 — > \fH and the dispersion relation contains only a single length scale 

(O; 

(5| - u 2 )ei {3 ( e-i - (bid 2 - u 2 - 1)} f h = 0. (76) 

It is possible to show analytically that the only root of the dispersion relation is Cj — 0, which means that dissolution 
is neutrally stable in the diffusive limit (H — > 00). On the other hand, the numerical results in Fig. [3] imply that the 
dissolution is unstable for even an infinitesimal convective flux. 



C. Geophysical implications 

Figure |4] summarizes the most important results of this study. Here we plot the (dimensionless) wavevector and 
growth rate of the dominant (most unstable) mode of the fracture instability. The convective limit extends up to 
H « 0.01; in this range both the dimensionless wavelength and growth rate are nearly constant. Thus for convection- 
dominated infiltration, the wavelength and timescale are simply related to the underlying geophysical parameters: 

\ ^^1* _ l-6fc e//7 

K e ff ho 

In order to put these results in a geological context, we consider typical values of the physical parameters charac- 
terizing dissolving fractures. Fracture apertures are between 0.005 cm and 0.1cm [2(J[28|,|2!|, and hydraulic gradients 
are of the order of 10~ 3 to 10 _1 [10, [3l|. This gives a range of characteristic flow velocities in undissolved frac- 
tures from 10 _4 cms _1 to 10cms~ 1 . The corresponding Peclet numbers are 10 _1 < Pe < 10 5 , taking the solute 
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FIG. 3: Growth rates of the inlet instability in the reaction-limited case (G = 0). The solid line corresponds to the purely 
convective case (H = 0). The other lines show results for increasing H: H = 0.1, H = 1, and H = 10. The dimensionless 
growth rate totd is plotted against the dimensionless wavevector 2-k/kX. 



diffusion coefficient D — 10 _5 cm 2 s _1 . The reaction rates vary widely, depending on the mineral. For example, 
relatively fast dissolving gypsum has a reaction rate k of the order of 0.01 cms -1 [32j, whereas siliceous minerals 
have surface reaction rates of the order of 10 -9 cms -1 [13, H3|. The typical reaction rates for calcite are in the range 
10 _5 cms _1 — 10 _4 cms~ 1 [13, Hl|. Thus the limitations imposed by the diffusion of reactant across the fracture 
aperture vary widely, resulting in a broad range of possible G values: from G ~ 10~ 7 in quartz, through G ~ 0.1 for 
a typical calcite fracture, up to G ~ 1 — 10 in gypsum. Nevertheless, Fig. [4] shows that both the maximal growth rate 
and the position of the maximally unstable wavelength depend only weakly on G; Cj max changes by 25% in the range 
< G < oo, with a similar change in the corresponding wavelength. However, the data in Fig. 0]is given in terms 
of dimensionless quan tities and the absolute growth rates vary dramatically across different minerals. For quartz, 
with 7 = 6- 10~ 5 [3(3], the time unit td ~ 5000 years, whereas the relevant timescale for calcite is a few months [34| . 
The same holds for the instability wavelengths, A, which vary from centimeters (gypsum) to kilometers (quartz). It 
is important to realize that the initial instability wavelength will in general be different from the spacing between 
protrusions in a mature formation. This is due to a coarsening of the pattern that is characteristic of this kind of 
dynamics 35]; the fingers compete with each other for the flow such that the longer ones grow more rapidly but the 
shorter ones become stagnant. As a result, the characteristic length between active (growing) protrusions increases 
with time. 

In geophysical systems, diffusion has only a small effect on the instability. Although H can vary from ~ 10 -15 (for 
wide fractures in siliceous formations) up to about 1 for narrow fractures in gypsum, fracture dissolution is typically 
convection dominated (H < 1). The residual diffusion leads to a slight shift of the peak growth rate towards longer 
wavelength, as observed in Fig. SJ but the wavelength and growth rate depend primarily on Da e // (j2"Tj) , via the 
penetration length l p and the dissolution time scale td, with just small corrections from H. 

These considerations refer to fracture dissolution in a natural geological setting. For carbonate acidization (e.g. 
with hydrochloric acid) the corresponding reaction rates are significantly higher than for dissolution with aqueous 
CO2; in acidization k ~ 10 _1 cms _1 [36], so that G can be larger than 100 (for ho ~ 0.1cm), which means that the 
dissolution rate is strongly limited by diffusion across the aperture. In the transport limit (G S> 1), H = Sh/Pe 2 is 
small under the typical flow rates used in acidization. 
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FIG. 4: Wavelength (Xmax) and growth rate (u} ma x) of the most unstable mode. The solid lines correspond to the exact solution 
of the reaction limit (1751) . The symbols indicate numerical results from the spectral code for a range of G: G = (circles), 
G = 1 (squares), and G — 00 (triangles). The dimensionless wavevector, u = 2n / kA, and growth rate, to = u>td, are plotted 
against the convective parameter H . The results for G = agree to at least 6 figures with the analytic dispersion relation from 
(|75|) . Up to 320 basis functions were needed to obtain convergent solutions (6 figures) from the spectral code at large H, but 
as few as 20 are sufficient for H < 1. 



D. Reaction order 

Experiments on the dissolution of limestone suggest that, near saturation, dissolution follows a nonlinear rate law, 
c.f. Eq. (Uni): 

R(c ca ) = -kc sat (l - — ) , n > 1, (78) 



c 



sat 



where c sa t is the saturation concentration of calcium ions. If we define a relative undersaturation c = (c sa t — c ca )/ (c sa j- 
Ci n ), where Ci n is the concentration of calcium ions at the inlet, then the transport equation, from is 



a. • 



q x d x c + q v d y c=-2k 1-— c n . (79) 

\ Csat J 

For simplicity, we only consider reaction-limited, convection-dominated dissolution. The equation describing aperture 
opening, analogous to (fTTj) . is 

\ 71-1 

Ci, ' 



d t h = 2kj 1 - — c", (80) 

where 7 = (c sat — Ci n ) / vc so i. The remaining equations, continuity and compatibility, are given by Eqs. (|26[) and (|27l) . 
Assuming the aperture in the base state is uniform, hb(x) = ho-, the base concentration profile is 

*M d i 2Hn-l)(l-c m /c sat r- 1 x \^ r _! 1 

Cftfacj = I H = l + fl — ri )«x 1 11 (81) 

V qo J 

where k — 2fcn(l — Ci„/c sa t) n_1 /(7o- In the limit n — > 1, Eq. (|8ip approaches the exponential base profile for linear 
reaction kinetics (fT5|) and the expression for k reduces to Eq. ([2~T|) . 

A dispersion equation for the growth rate can be obtained for non-linear kinetics by following the procedure in 
Sec. IIV1 starting with the analogues of Eqs. (I34|) - (l35l) : 

Sq x d x c b + qbd x Sc = --d t Sh, (82) 

7 



V C S at J 



71-1 

71-1 



cW/i = 2fc 7 1 - — nc""Mc. (83) 
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FIG. 5: The impact of kinetic order on the growth rates. The growth rate of the instability is shown for various powers of n, 
including the linear rate law, n — 1, and the limit of high reaction order, n — > oo. 



The continuity and compatibility relations are the same as Eqs. (|36j) - (|37[) . Introducing dimensionless variables: 

£ = kx, V = Ky, t — , (84) 

ho 

and scaling Sh and q as in (|3T)]) . we obtain the following equations for f c , fh, and f q , denned in Eqs. (l4"T)) - (|4l?)) : 

f q dtSb + dtf c = -~f h , (85) 



/J 



-h = erVc (86) 



n 



(9| - u 2 )f q = -3u 2 f h . (87) 

The inlet saturation, Cj n , has been absorbed into the length and time scales (|84p. 

The base concentration (£&) can be eliminated from the equations for transport (1851) and erosion (|86|) by using (|8ip : 



/, = [1 + (1 - n" 1 )?] (<&/a + n^/ c ) 
cj// l = n[l + (l-n- 1 X] _1 /c- (89) 
Combining these equations with (1871) we get a dispersion equation for arbitrary kinetic order, 

n 

u>{8$ - u 2 ) [1 + (1 - n~ l )(\ »-i [1 + 9 5 (1 + (1 - n- l )i)) f h + 3u 2 f h = 0, (90) 

which is well behaved in the limits n — > 1 and n — ¥ oo. 

The impact of kinetic order is illustrated in Fig. [SJ which shows that even strongly non- linear reaction kinetics (n — > 
oo) do not suppress the instability. The dimensionless growth rate depends only weakly on reaction order, reflecting 
our choice of scaling for the dimensionless length and time. Thus, as a first approximation we can take the peak growth 
rate as Cj max ~ 1 and the corresponding wavevector u max ~ 1, independent of reaction order. Then, in absolute terms, 

the wavelength corresponding to maximum growth is roughly proportional to Xmix ~ A ma3; /n(l — Ci n / Csat) 11 " 1 , 
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FIG. 6: The effect of finite system size on the growth rate in the convection-reaction (G = H = 0) limit. 



where X ma x ~ 27rgo/2fc is the peak wavelength for linear kinetics. This is slightly counterintuitive since increasing 
reaction order tends to increase the penetration of reactant into the fracture. Nevertheless its effect on the instability 
is to shorten the wavelength of the most unstable mode. However the wavelength is also strongly dependent on Cj n , 
and a partially saturated solution at the inlet increases the wavelength of the most unstable mode. The inlet solution 
to the fracture must be nearly saturated (cj„ — > c sat ) for non-linear kinetics to apply 18], so the wavelength in such 
cases is almost entirely dependent on the extent of the (small) undersaturation. The corresponding growth rate of 

the instability u« = w TO ax(l — Cinl Csat) n ~ x is sharply limited by the degree of undersaturation. 



E. Finite length fractures 



The previous analysis corresponds to a semi-infinite system, x > 0, which is the relevant limit for geophysical 
systems where the length of the system, L, is usually many orders of magnitude larger than the penetration length . 
However, in laboratory experiments as well as in petroleum reservoir stimulation, the relevant length scales are much 
smaller and finite-size effects may be important. In this case, the far-field boundary condition q x (x — ¥ oo, y,t) = qo 
must be replaced by a constant pressure condition at the outlet; then q y (x = L, y, t) = or, in terms of perturbations, 

Sq v (x = L) = 0. (91) 

Figure [6] shows the effect of a finite length aperture in reaction limited, convection-dominated dissolution (H = G = 
0). Now all three solutions from Eq. (|73p are needed; Eqs. (|5cj|) and (|i?Tj) fix the perturbation to within an arbitrary 
amplitude, while Eq. (|59j) enforces the eigenvalue condition. The additional length scale leads to a richer spectrum of 
possibilities; in particular, the longest wavelengths are now less stable than in unbounded [L — > oo) fractures. The 
shape of the dispersion curve changes considerably as the length of the system is reduced and for short fractures, 
(kL < 2), the growth rate is maximum at zero wavevector. As the length of the fracture increases, the wavelength of 
the most unstable mode shifts to larger u and the longest wavelengths are only weakly unstable; as L — > oo the growth 
rate at zero wavevector vanishes altogether. In fact, the growth rate at u = has a particularly simple analytical 
form 

3 [l - (1 + nL)e- KL ] , 
u>(u - 0) = -i ^ ' i, (92) 
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FIG. 7: Growth rate in the long-wavelength (u — > 0) limit, loo, in the convection-reaction (G — H = 0) limit. 

which is shown in Fig. [Jj Both for very small and very large lengths the long-wavelength growth rate is relatively 
small, with a maximum at kL sa 1.8. 

An analysis of Fig. [6j together with Fig. [2] offers some insight into the typical dispersion curve for a fracture 
dissolution instability, which exhibits a strong wavelength selection with a well-defined maximum in the growth rate 
for X max ~ kT 1 . The results presented in this section show that stabilization of the growth of long wavelength 
instabilities is connected with the far-field boundary condition (|60p. which imposes a uniform flow at large distances 
from the inlet. However, in a finite system, the constant pressure condition at x — L does not require q x to be 
uniform, and hence does not lead to a stabilization of long- wavelength modes, as shown in Figs. [5] and [7J On the other 
hand, Fig. [2] shows that the shape of the short- wavelength spectrum is controlled by reaction kinetics. In particular, 
transport-limited kinetics decreases the short-wavelength growth rates, since in this regime dissolution slows down as 
the fracture opens (J8]). 

VII. CONCLUSIONS 

In this paper, we have analyzed the stability of a one-dimensional reaction front in dissolving fractures. Strikingly, 
the dissolution front turns out to be unstable over a wide range of wavelengths, suggesting that fracture dissolution 
is an inherently two-dimensional process. The maximal growth rate corresponds to wavelengths of the order of the 
penetration length n~ l and this result turns out to be remarkably insensitive to the details of the reaction and 
transport mechanisms in the fracture: the maximum is shifted towards longer wavelengths when strong diffusion is 
present or for strongly nonlinear reaction kinetics, but the shift is relatively small and n\ max remains within the same 
order of magnitude. The only case where there is a qualitative change in the dispersion curve is a finite-length system. 
For relatively short fractures, kL < 3, the maximum growth rate occurs at zero wavevector and long-wavelength 
modes remain unstable. 

In summary, the reactive front instability has been shown to be a generic phenomenon in the dissolution of fractured 
rock. Hence the predictions of fracture breakthrough times, crucial for speleogenesis and for the assessment of 
subsidence hazards, cannot be based on one-dimensional models. Instead, a two-dimensional model is necessary to 
take into account the highly localized dissolution front. Numerical [H, 0, S] and theoretical [35[ work has suggested 
that the dissolutional instability leads to a strong focusing of the fluid flow into a few active channels, which advance 
in the fracture while competing with each other for the available reactant. However, a quantitative characterization 
of this non-linear process, which is essential for the prediction of fracture breakthrough times, remains elusive. 
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Appendix A: Convection-diffusion equation for reactant transport 

In this Appendix we will investigate the validity of the two-dimensional steady-state transport equation J3J in 
various parameter ranges. In the spirit of the Reynolds approximation, we will then assume that the global solution 
for parallel plates can be applied locally, if the fracture aperture field varies sufficiently slowly, |V/i| ~ 0{1). Taking 
the flow to be along the x-axis and the normal to the fracture surfaces along the z-axis, the convection-diffusion 
equation for the three-dimensional concentration field c 3 d{x, z,t) can be written as 

d t c 3d + v x d x c 3 d = D (d 2 c 3d + d 2 c 3d ) , (Al) 

where v x = 6v a (z/h — z 2 /h 2 ) and v a = q x /h is the aperture-averaged fluid velocity. In addition we have boundary 
conditions on the fracture surfaces 

Dd z c 3d \ z =o = kc 3d , Dd z c 3d \z=h = -kc 3d , (A2) 

and at the inlet, 

C3dU=0 = Qti, (A3) 

where Cj„ is the inlet concentration. 

A direct integration of (|A1|) over the z coordinate gives a two-dimensional averaged convection-diffusion-reaction 
equation involving three different concentrations, 

2k 

d t c a + v a d x c = Dd 2 c a - — c w ; (A4) 

h 

c a (x,t) = hr 1 C3d(x, z,t)dz, is the aperture-averaged concentration, c is the cup-mixing concentration ([3]), and 
c w (x,t) = c 3 d{x,0,t) — C3d(x,h,t) is the reactant concentration at the fracture surfaces. Following standard proce- 
dures for averaging the convection-diffusion equation, we will solve Eqs. (|A1|) - (IA3|) to find relations between these 
average concentration fields in different parameter ranges. In particular we will show that Eq. ([3]) is correct in the 
important limits of convection-dominated transport (Sec. IA 2[) and reaction-limited (Sec. IA 3[) kinetics. 



1. Scaling and steady state 

The steady state approximation in ((3]) can be justified by the time-scale separation between the transport of reactants 
and the consequent change in fracture aperture. The dissolution time scale is characterized by id — h/2k-f — tdj (1 + G) 
(|12p . where the acid capacity number 7 = Ci n /vc so i is usually small, because of the high molar concentration of the 
solid phase. For example calcite contains roughly 25 moles per liter, whereas even a strong acid is rarely used in more 
than 1 molar concentrations; in the natural dissolution of calcite by atmospheric CO2, 7 ~ 10 -4 . To see how a small 
7 leads to the steady-state limit we scale the time by Id in addition to the usual scaling of lengths: 

«-ff-S' T -c- <A5) 

The axial distance is scaled by the characteristic length I = v a h/2k, and the transverse distance is scaled by h. In 
addition the fluid velocity is scaled by v a and the concentration by Ci n : 

£ C = - = 6C-6C 2 , c 3d =°^. (A6) 

The scaled convection-diffusion equation, 

jd T c 3d + V£d 6 c 3d = Hd^c 3d + G^dl^d, (A7) 
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is then be characterized by 7 and two new dimensionless groups: G = 2kh/D and H — 2kD/v 2 h. G and H are related 
to the corresponding parameters defined in the main body of the paper by G = GSh ([13]). and H = H(l + G) (fP51) . 
In this appendix we consider the transverse (z) direction explicitly and so the Sherwood number does not appear 
in the defining equations; the ratio of diffusive and reactive fluxes is then characterized by G rather than G. Since 
the reactive flux appears in the boundary conditions rather than the underlying equations, the ratio of diffusive and 
convective fluxes is more naturally defined by H rather than H. 
The steady-state convection-diffusion equation 

v&c 3d = Hd\hd + G^dl^d (A8) 

is reached in the limit 7 — > 0, and is valid under most circumstances arising in fracture dissolution. The boundary 
conditions in the dimensionless variables are 

d ( c 3d = ±^c 3d , C = 0,1; (A9) 

and the average equation for steady-state reactant transport is 

d^c = H3lc a -c w . (A10) 



2. Convective limit: H = 0. 



Fracture dissolution is usually characterized by small H, corresponding to the convective limit H — > (Sec. IVI Cl) . 
The time- independent convection-diffusion equation is then 



v&c 3d = G~ L d z c c 3d , (All) 

which can be solved by separation of variables, c 3 d = /(05(C) [HI- The decay in the axial direction is a sum of 
exponentials, exp(— A„£), where A„ are related to the positive eigenvalues of the equation 

d 2 c g + 16r 2 ((-( 2 )g = 0, (A12) 



with r n = y 3GA„/8. This equation has a single solution that satisfies the symmetry condition g(0) = g(l), 

9(0 = if! i; r(2( - l) 2 ) e" 2 ^ 1 ). (A13) 

Applying the boundary conditions from (|A9|) leads to the eigenvalue equation for r(G), 

r(r-l) 1*1 + ( r_ f) lFl (^T^ ;r ) (A14) 

The average equation for the concentration, 

<%c = ~c w , (A15) 

implies that for a single mode Ac = c w (the same result follows from integrating Eq. (|A12|) over Q. Using the Sherwood 
number to connect c w and c ((6]), 

(A16) 



1 + G/Sh' 

we can relate the eigenvalue A = 8r 2 /3G to Sh 



Sh = ^. (A17) 



Thus the steady-state convection-reaction equation is simply 



18 



where Sh(G) is determined from the smallest root of (IA14[) . 

For reaction-limited kinetics r — > 0, and the hypergeometric functions in Eq. (|A14[) can be expanded around r = 0; 
solving for G we obtain a quadratic equation for A, 

G = \r 2 + gr 4 + 0(r 6 ) = AG + i^A 2 G 2 , (A19) 

with a solution A = 1 — 17G/140 4- C(G 2 ). The concentration is nearly uniform across the aperture and decays 
axially as a single exponential e~ A ^. From Eq. (|A17[) we find the Sherwood number for reaction-limited kinetics 
Sh° = 140/17 « 8.24. 

In the transport limit the concentration at the walls vanishes (Graetz problem) and the eigenvalues A„ = 8r 2 /3G 
can be found from the roots of the equation 

iFi (V"^ ;r ) =o - (A2o) 

The transport-limited Sherwood number, Sh°° w 7.541, follows from the smallest eigenvalue r ~ 1.6816. In the 
numerical work we will ignore the weak dependence of Sherwood number on G and take Sh = 8 throughout. 



3. Reaction-limit: G — > 0. 



Away from the convective limit, the diffusive flux prevents a solution of the transport equation (|A8|) by separation 
of variables. However, when the reaction rate is small, such that G <C 1, the deviation in concentration from the 
average concentration, c^d — c Q , can be expanded in powers of G [37l l38j. 

c 3 d -Ca = Gc« + G 2 c< 2 ) + . . . ; (A21) 

it follows that 

l 

c {l) d( = 0. (A22) 

From Eq. (|A8|) , the zeroth order convection-diffusion equation is 

(6C - 6C 2 )<%c a = Hd^ca + <9 2 c (1) . (A23) 

Integrating Eq. (|A23|) across the aperture and using the boundary condition (|A9|) d^c^ = ±c a /2, we obtain the 
average equation 

8 s c a = Hd}c a - c a , (A24) 

which is the reaction limit of (|A10I) . In this limit the concentration profile is uniform across the aperture and all three 
concentrations, c, c a , and c w are equal. 

Equation (IA24p can be subtracted from (|A23[) to eliminate the diffusion term, 

(6C - 6C 2 - l)%c -ca = d 2 c c^ . (A25) 

Solving for 

cm = e(i J e -l?-l? + ^-tJ£-i + l-)., (A26) 



the linear term in ( is introduced to satisfy the boundary conditions in (|A9[) and the constant term is to enforce the 
condition in (|A22I) . Finally, we use Eq. (|A26[) to relate c and c w to c Q : 



C = J\&C - K 2 Wa + GC«R = f 1 + ^ j C a - ^0 ( C a , 



(A27) 



c a + Gc«(C = 0) = [ 1 - ^ ) c a + ^c a . (A28) 
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These are the equivalents of the results in [38| (67c & d), but for flat plates instead of tubes. 

Using Eqs. (|A27[) and (|A28p to eliminate c a and c w from the average equation (|A10[) . the transport equation 
becomes 



The third-order term in Eq. (|A29|) . 



is small compared with the convective term, 



HG „, v a h h 2 „„ , . s 

dfc = die, A30 

210 ? 2k 210 x K ' 



d 6 c = ^-d x c, (A31) 

on all scales larger than the aperture h. Since h is small on scales of interest in fracture dissolution we can safely 
ignore this term. Similarly, the diffusive term (4#G/105)<9|c is small compared to c. Dropping these terms leaves the 

renomalization of the reaction term as the leading-order correction for finite G (in the steady-state limit), 

d € c = Hdic 1 jr. (A32) 

5 C 1 + G/Sh° V ; 

The average equation for the cup-mixing concentration has no Taylor dispersion term, but only the contribution from 
molecular diffusion. This is true both in the convective limit (arbitrary G) and the reaction limit (arbitrary H). 



4. Summary 

In this appendix we have examined the structure of the depth-averaged convection-diffusion equation across a 
range of Damkohlcr and Peclet numbers. The dimensionless parameter H = Da e ///Pe is usually small in fracture 
dissolution, which implies a convection-dominated process. In such cases the steady-state convection-reaction equation 
(|A18[) follows (see Sec. IA 21) . with only a weak dependence of the Sherwood number on reaction rate and entrance 
length. 

When diffusion plays a significant role, the structure of the average equations is more complex, and it is not 
possible to rigorously treat transport in the case of sig nificant transverse and axial diffusion (G ^> 1, H ^> 1) without 
considering more than one average concentration [37l l38j. Nevertheless, in Sec. IA 3l we showed that in the reaction 
limit (G < 1) the structure of Eq. (|A18|) is preserved (|A32|) . 



Appendix B: Scale-dependent Peclet and Damkohler numbers 

The one-dimensional transport equation ()13[) can be non-dimensionalized by the penetration length 

q oK d ( c - Dh oK 2 dlc = (Bl) 
where £ = kx. Dividing Eq. (|B1[) by qgK suggests two new dimensionless constants: 

^ eK ~ n u - ~T~ ' D&K — /-, ; n \ — — u — ' V B2 J 

Dung nriQ 9o K (l + G) nlio 

Pe K is the ratio of convective to diffusive fluxes on the length scale while Da K is the ratio of convective to reactive 
fluxes on the same scale; Da K is based on the effective reaction rate k e ff ©. The transport equation on the scale of 
the penetration length k _1 is then 



d^c — Pe K 1 9|c = — Da K c. 



(B3) 
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The parameter H retains the same meaning with the new definitions of Peclet and Damkohler number, 

Da^Da^ 

Pc Pe K ' V ; 

and the two new parameters can be written solely in terms of H: 

2 2H 
Pc K = , , Da K = , . (B5) 

Although Pe K and Da K are not independent, Da K = 1 -l-Pe" 1 , it is a notational convenience to treat them so; however 
the results are discussed in terms of the independent parameters G and H . 

On the relevant length scale for fracture dissolution, the ratio of convective and diffusive fluxes is characterized 
by Pe K . Nevertheless we prefer to characterize the dissolution in terms of G and H rather than G and Pe K , since 
both Pe K and Da K have simple expressions in terms of H. In the convective limit (the most important for fracture 
dissolution) H — > Pe" 1 , while in the diffusive limit H — > Pe~ 2 . Thus the convective limit implies Pe K — > oo and 
H — > 0, while the diffusive limit is the opposite, but the mapping is not a simple inverse relation. 



Appendix C: Derivation of the compatibility relation 

Throughout the paper we will frequently make use of the compatibility relation (|27[) . which can be derived by noting 
that, from fTJ), 

d y q x = 

Similarly 

d x q y = 

Subtracting (|C2[) from (|C1[) leads to the compatibility relation 

3 3 
d y q x ~ Tlxdyh = d x q y - -q y d x h. (C3) 



h 3 d x „p — 3 h 2 d„hd x p 

Yl\x 12^ 



1 3 
— h 3 d xy p+ -q x d y h. 



(CI) 



h 3 d xv p — 3 h 2 d x hd v p 



1 3 
—h d xy p + -q y d x h. 



(C2) 
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